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Abstract 

For  Stokes  equations  with  discontinuous  viscosity  across  an  arbitrary  interface  or/and 
singular  forces  along  the  interface,  the  pressure  is  known  to  be  discontinuous  and  the 
velocity  is  known  to  be  non-smooth.  It  has  been  shown  that  these  discontinuities  are 
coupled  together  which  makes  it  difficult  to  obtain  accurate  numerical  solutions.  In  this 
paper,  a  second  order  accurate  numerical  method  that  decouples  the  jump  conditions  of 
the  fluid  variables  through  two  augmented  variables  has  been  developed.  The  GMRES  it¬ 
erative  method  is  used  to  solve  the  Schur  complement  system  for  the  augmented  variables 
which  are  only  defined  on  the  interface.  The  augmented  approach  also  rescales  the  Stokes 
equations  in  such  a  way  that  fast  Poisson  solvers  can  be  used  in  each  iteration.  Numeri¬ 
cal  examples  against  exact  solutions  show  that  the  new  method  has  average  second  order 
accuracy  in  the  infinity  norm,  and  the  number  of  GMRES  iterations  is  independent  of 
mesh  sizes.  An  example  of  a  moving  interface  problem  is  also  presented. 

Keywords:  Incompressible  Stokes  flow,  discontinuous  viscosity,  coupled  jump  conditions,  in¬ 
terface,  discontinuous  and  non-smooth  solution,  immersed  interface  method,  fast  Pois¬ 
son  solver 
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1  Introduction 

The  incompressible  Stokes  or  Navier-Stokes  equations  with  discontinuous  viscosity  and  sin¬ 
gular  forces  arise  from  many  important  applications  in  fluid  and  biofluid  mechanics.  One  of 
particular  example  is  Peskin’s  immersed  boundary  (IB)  model  to  simulate  the  blood  flow  in  a 
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human’s  heart  [19,  20,  23].  The  idea  of  IB  formulation  is  to  treat  the  complicated  immersed 
boundary  (such  as  a  heart  valve)  as  a  force  generator  in  the  fluid  domain,  or  mathematically, 
a  Dirac  delta  function  force  distribution  along  the  immersed  boundary.  While  Peskin’s  IB 
model  has  been  intensively  studied  in  the  literature,  almost  none  has  studied  the  case  where 
the  viscosity  is  discontinuous.  In  this  paper,  we  consider  the  following  two-dimensional  sta¬ 
tionary  incompressible  Stokes  equations 

Vp  =  V  •  n  (Vu  +  (Vu)T)  +  g  +  '  f(s)  ^(x  -  X(s))ds,  x  £  Cl,  (1.1) 

V  •  u  =  0,  x£(l,  (1-2) 

where  u  =  (tt,  v )  is  the  velocity  vector,  p  is  the  pressure,  T  is  an  arbitrary  interface  parame¬ 
terized  by  the  arc-length  s ,  f  is  the  density  of  the  force  strength  along  T1,  p  is  the  viscosity, 
which  we  assume  to  be  a  piecewise  constant 


h  = 


if  x  £  D+, 
if  x  £  D~, 


(1.3) 


with  a  finite  jump  across  T,  g(x)  is  a  bounded  forcing  term,  e.g.,  the  gravitational  force, 
which  can  be  discontinuous  across  T  as  well,  and  D  is  a  bounded  domain  which  we  assume 
to  be  a  rectangle,  see  Fig.  1  (a)  for  an  illustration  of  the  problem.  We  assume  periodic 
boundary  conditions  for  the  pressure  and  the  velocity  in  this  paper.  The  system  (1.1)-(1.3)  is 
the  model  equations  that  we  are  going  to  discuss  in  this  paper.  The  existence  of  the  solution 
to  the  system  (1.1)-(1.2)  can  be  found  in  [22]. 

There  are  at  least  two  difficulties  in  solving  (1.1)-(1.2)  numerically.  The  first  one  is  to  deal 
with  the  singular  force  term.  A  simple  way  is  to  use  Peskin’s  discrete  delta  function  approach 
to  distribute  the  singular  force  to  nearby  grid  points.  Such  a  discretization  is  typically  first 
order  accurate  and  will  smooth  out  the  solution.  The  second  difficulty  is  how  to  deal  with 
the  discontinuous  viscosity.  A  simple  smoothing  method  may  have  large  errors,  see  [15]  for  a 
one  dimensional  example  there. 

In  the  case  of  the  continuous  viscosity,  various  methods  have  been  developed  in  the  litera¬ 
ture.  We  refer  the  readers  to  [3,  4,  12,  18,  23]  for  various  methods  and  the  references  therein. 
The  difficulty  with  discontinuous  viscosity  is  that  the  jump  conditions  for  the  pressure  and 
velocity  are  coupled  together,  see  (2.4)-(2.7),  which  makes  it  difficult  to  discretize  the  system 
accurately. 

In  this  paper,  we  develop  a  new  second-order  sharp  interface  method  which  uses  the  exact 
jump  conditions  for  two-phase  incompressible  Stokes  flow  with  discontinuous  viscosity.  We 
believe  that  our  method  is  the  first  sharp  interface  method  with  second-order  accuracy.  The 


1The  singular  source  term  fr  f(s)  ^(x  —  X(s))ds  can  also  be  written  as  ((f  ■  n)n  +  (f  ■  r)r)  (5(F),  or  in 
the  form  of  ((f  •  n)n  +  (f  ■  t)t)  |Vy>|  where  ip  is  a  level  set  function  whose  zero  level  set  represents  the 
interface  T. 
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(a) 


(b) 


A 
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Figure  1:  (a):  A  diagram  for  the  incompressible  Stokes  equations  defined  on  a  domain  17  with 
an  interface  T  across  which  the  viscosity  //  is  discontinuous,  (b):  Force  density  decomposition 
in  which  f\  and  /2  are  the  force  density  in  the  x-  and  y-  directions,  while  /i  and  f-2  are  the 
force  density  in  the  normal  and  tangential  directions. 

idea  is  to  introduce  two  augmented  variables  that  are  defined  only  along  the  interface  so  that 
the  jump  conditions  can  be  decoupled  and  the  immersed  interface  method  [11,  13]  can  be 
applied.  The  GMRES  iterative  method  is  applied  to  the  Schur  complement  system  for  the 
discrete  augmented  variables.  Furthermore,  our  approach  rescales  the  original  problem  and 
enables  us  to  use  fast  Poisson  solvers  in  the  iterative  process.  Each  GMRES  iteration  requires 
to  solve  the  rescaled  Stokes  equations,  which  can  be  done  by  calling  a  fast  Poisson  solver  three 
times,  and  an  interpolation  scheme  to  evaluate  the  residual  of  the  Schur  complement. 

We  should  mention  other  related  work  in  this  area.  One  is  the  fast  Stokes  solver  on 
irregular  domains  based  on  an  integral  equation  approach  [1,  5]  which  may  be  applied  to  the 
two-phase  Stokes  flow  discussed  here.  There  are  also  some  work  on  the  Navier-Stokes  equa¬ 
tions  with  discontinuous  viscosity,  for  example,  the  first  order  accurate  ghost  fluid  method 
[9].  The  system  of  Navier-Stokes  equations  is  a  time  dependent  problem  in  which  we  can  use 
a  time  marching  method.  The  Stokes  flow  discussed  here,  on  the  other  hand,  is  an  elliptic 
system  in  which  we  have  to  solve  the  entire  system  simultaneously. 

The  paper  is  organized  as  follows.  In  the  next  section,  we  give  some  discussions  on  the 
jump  conditions  of  the  Stokes  equations  involving  an  interface.  We  will  see  how  the  jump 
conditions  for  the  pressure  and  velocity  are  coupled  together,  and  explain  the  idea  how  to 
decouple  the  jump  conditions  by  introducing  two  augmented  variables  and  two  augmented 
equations.  In  Sec.  3,  we  present  the  new  algorithm  in  detail.  We  explain  how  the  GMRES 
iterative  method  can  be  applied  to  the  discrete  augmented  variables  without  explicitly  form 
the  coefficient  matrix.  In  Sec.  4,  we  present  some  numerical  experiments  against  the  exact 
solutions  to  check  the  performance  of  new  method.  An  example  of  moving  interface  is  also 
presented  there. 
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2  Jump  conditions  for  incompressible  two-phase  Stokes  equa¬ 
tions  with  discontinuous  viscosity  and  singular  forces 

Referring  to  Fig  1  (a)  and  the  equations  ( 1. 1)-(  1 .3) ,  we  assume  the  interface  T  is  a  smooth 
curve.  At  a  point  (X,  Y )  on  the  interface,  we  use  the  notations  n  and  r  to  represent  the  unit 
normal  and  tangential  directions  respectively.  We  use  /i(s)  =  f(s)  •  n  and  f2(s)  =  f(s)  •  r 
to  represent  the  force  density  in  the  normal  and  tangential  directions.  If  T  is  smooth, 
then  in  a  neighborhood  of  T,  the  distance  function  d(x,T)  is  also  a  smooth  function.  The 
normal  and  tangential  directions  of  T  can  be  extended  to  the  neighborhood,  for  example 
n  =  Vd(x,  r)/|Vd(x,  r)|.  Therefore  the  quantities  such  as  u  n,  etc.,  are  well  defined 

in  the  neighborhood  of  T. 

In  [8] ,  the  interface  conditions  for  two-phase  Stokes  equations  with  discontinuous  viscosity 
and  a  Dirac  delta  function  singularity  are  derived  which  is  summarized  in  the  following 
theorem: 


Theorem  1  Assume  T(s)  £  C2,  fi(s)  £  C1,  and  f-2(s)  £  C1.  Let  p  and  u  =  (u,v)  be  the 
solution  to  the  Stokes  equations  (1.1)-(1.2).  We  have  the  following  jump  conditions  across 
the  interface  T. 


\p\  =  2 

dp 
dn 

du 


<9u 

hw-  ■  n 
dn 


fi, 


r  ,  3  i 

d2 

=  [g  •  n]  +  — /2  +  2 

h~Q^2  (u'n) 

+ 


du 


+  h  =  o, 


(2.4) 

(2.5) 

(2.6) 


[//V  •  u]  =  0. 


(2.7) 


where  the  jump  [  •  ]  of  a  quantity,  for  example,  [p]  at  a  point  X  is  defined  as 


\P] 


lim  p(x)  —  lim  p(x). 
xer  x-+x.xefi+  x->x.xeo- 


(2.8) 


We  will  omit  the  subscript  X  £  T  in  the  rest  of  the  paper  for  simplicity  of  the  notations. 
The  first  jump  condition  (2.4)  is  the  result  of  balancing  force  in  the  normal  direction  while 
the  third  one  (2.6)  is  the  result  of  balancing  force  in  the  tangential  direction.  The  last  jump 
condition  (2.7)  is  obtained  from  the  incompressibility  condition.  The  jump  condition  for  the 
normal  derivative  of  the  pressure  (2.5)  can  be  obtained  by  applying  the  divergence  operator 
to  the  momentum  equation  (1.1),  see  [8].  It  is  worthy  mentioning  that  the  equation  (1.1)- 
(1.2)  can  be  re-written  as  the  one  without  the  singular  force  term  but  accompanied  with  the 
above  jump  conditions. 
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From  the  incompressibility  condition,  one  can  easily  prove  that  •  n]  =  0.  If  p  is 
continuous,  then  •  n]  =  0  and  [p  ( u  •  n)]  =  0,  and  the  jump  conditions  for  the  pressure 
and  the  velocity  are  decoupled  from  (2.4)-(2.7).  In  this  case,  we  recover  the  jump  conditions 
derived  and  used  in  [12,  10].  The  second  order  accurate  immersed  interface  method  has  been 
developed  in  [12,  13]  if  the  viscosity  is  continuous.  When  p  has  a  finite  jump  across  the 
interface  T,  to  our  best  knowledge,  there  is  not  yet  second  order  accurate  sharp  interface 
method  in  the  literature  for  the  Stokes  equations.  To  get  a  second  order  accurate  algorithm 
based  on  the  immersed  interface  method,  our  strategy  is  to  introduce  two  intermediate, 
or  augmented  variables,  along  the  interface  so  that  the  jump  conditions  can  be  decoupled. 
Meanwhile,  we  also  need  two  augmented  equations  to  complete  the  system  of  the  governing 
equations. 

There  are  more  than  one  ways  to  introduce  augmented  variables  so  that  the  jump  con¬ 
ditions  can  be  decoupled.  Different  augmented  variables  and  equations  will  lead  to  different 
algorithms.  Note  that  for  discontinuous  viscosity,  there  are  at  least  two  different  scales  corre¬ 
sponding  to  the  two-phase  flow.  Furthermore,  with  constant  viscosity,  we  can  use  fast  Poisson 
solvers  to  get  the  solution  of  the  Stokes  flow  [12,  13].  Based  on  these  two  considerations,  we 
introduce  the  jumps  [/xu](s)  and  [/xu](s)  along  the  interface  as  two  augmented  variables. 
The  advantages  and  details  can  be  seen  in  the  rest  of  the  paper. 

Before  we  proceed,  we  introduce  the  local  coordinates  transform  at  a  point  ( X ,  Y )  on  the 
interface  F  as, 

f  £  =  (x  —  X)  cos  9  +  (y  —  Y)  sin  9, 

\  (2.9) 

[  77  =  — (x  —  X)  sin  6  +  (y  —  Y)  cos  9, 

where  9  is  the  angle  between  the  x-axis  and  the  normal  direction  at  the  point  (X,Y),  see 
Fig.  1  (b).  Using  this  local  coordinate  system,  we  can  rewrite  the  last  two  jump  conditions 
(2.6)-(2.7)  in  terms  of  the  augmented  variables  [/xu](s)  and  [/iw](s)  as  follows. 


Lemma  1  Let  p,  u,  and  v  be  the  solution  to  the  Stokes  equation  (1.1)-(1.2).  Define 


U  =  pU,  V  =  pV,  u  =  (ft,  v). 


Then  we  have  the  following  jump  relations  for  u  and  v 


du 

dn 


n 


sin  9  — 


cos  9, 


dv 

dn 


du 

du 

if'2  + 

Jh  '  n 

I  cos  9  — 

Jh-  '  T 

du 

dn 


(2.10) 


(2.11) 

(2.12) 

(2.13) 


Proof:  Note  that  n  =  (cos  9,  sin  9)  and  r  =  (—  sin  9,  cos  9).  Re-write  the  incompress- 
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ibility  condition  [pA7  •  u]  =  0  in  the  local  coordinates,  we  have 


which  is 


du 

dn 


cos  9  — 


du 


sin  9  + 


dv 

dn 


sin  9  + 


dv 

frr 


cos  9  =  0, 


du 

dn 


cos  9  + 


dv 

dn 


sin  9 


Re-write  the  interface  relation  (2.6)  in  the  local  coordinates,  we  have 


du 

dv 

du 

dn 

sin  9  + 

dn 

COS  9  =  — f‘2  — 

— —  ■  n 

dr 

(2.14) 


(2.15) 


(2.16) 


From  the  two  equalities  (2.15)  and  (2.16)  above,  we  solve  for  [|^]  and  [|^]  to  get  (2.11)  and 
(2.12).  The  last  equality  is  verified  by  substituting  [|^]  with  (2.11),  and  [^]  with  (2.12)  in 
the  following 


<9u 

dn 


du 

dn 


cos  9  + 


dv 

dn 


sin  9 


3  The  numerical  algorithm 

Our  numerical  method  is  based  on  the  following  theorem. 

Theorem  2  Let  p,  u,  and  v  be  the  solution  to  the  Stokes  equations  (1.1)-(1.2).  Let  qi(s)  = 
[ft](s)  =  [pu](s),  q2(s )  =  [6](s)  =  \pv] (s),  and  q(s)  =  (qi(s) ,  q2(s)) .  Then  u,  v,  p,  q^s), 
q2(s)  are  the  solution  of  the  following  augmented  system  of  partial  differential  equations: 

(3.17) 

(3.18) 

(3.19) 


u 

h 


Ap  =  V  •  g, 

\p]  =  h~  •  t, 


dp 

dn 


9f2  092(q-n) 

dr  dr2  ’ 


Au  =  px  -  gi, 

[y\  =  qi, 

Av  =Py-  02, 

[®]  =  92, 


du 

dn 


dv 

dn 


f  3q  \  •  a  ( 5q  I  a 
f2  +  —  •  n  I  sin  9  -  1  —  •  r  )  cos  9, 


i  ,  <9q  \  ( dq 

/2  +  —  •  n  I  cos  0  -  I  —  •  r  ]  sm  9, 


(3.20) 
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The  proof  of  the  theorem  is  straightforward  from  the  Stokes  equations  (1. 1)-(1 .2) ,  the  jump 
conditions  in  Theorem  1  and  Lemma  1.  The  periodic  boundary  condition  is  used  so  that  we 
are  not  introducing  additional  boundary  condition  for  the  pressure. 

The  existence  and  uniqueness  of  the  solution  to  the  system  above  is  the  same  as  the 
original  incompressible  Stokes  equations  (1.1)-(1.2).  This  is  because  if  (u,v)  and  p  are  the 
solution  to  the  original  Stokes  equations,  then  they  are  also  the  solution  to  the  system  above 
according  to  the  definition  of  (u,v),  Theorem  1,  and  Lemma  1.  On  the  other  hand,  if  (u,v) 
and  p  are  the  solution  to  the  system  (3.17)-(3.20)  above  plus  the  periodic  boundary  condition, 
then  they  satisfy  all  the  equations  in  (1.1)-(1.2)  and  the  incompressibility  condition.  So  they 
are  also  solution  to  the  original  problem. 

Notice  that  if  we  know  q,  then  the  jump  conditions  for  the  pressure  are  all  known  and  we 
can  solve  the  pressure  independent  of  the  velocity.  After  the  pressure  is  solved,  we  can  solve 
the  velocity  from  (3.18)  and  (3.19).  The  three  equations  with  the  given  jump  conditions  can 
be  solved  using  the  immersed  interface  method  [11,  12,  13]  in  which  a  fast  solver  is  called 
with  modified  right  hand  sides  at  grid  points  near  or  on  the  interface.  This  observation  is 
the  basis  of  our  new  method.  We  can  try  to  find  q  iteratively  starting  from  an  initial  guess. 
The  compatibility  condition  for  the  two  augmented  variables  are  the  two  equations  in  (3.20) 
(which  means  that  the  velocity  is  continuous  across  the  interface).  It  is  also  important  to 
mention  that  the  incompressibility  condition  is  used  to  obtain  the  pressure  Poisson  equation 
of  (3.17). 

Once  the  augmented  variables  ([//it]  and  [pv])  and  the  augmented  equations  (the  two 
equations  in  (3.20))  are  chosen,  the  success  of  the  numerical  algorithm  depends  on  how 
efficiently  we  can  solve  the  augmented  variables.  Note  that  the  augmented  approaches  have 
been  developed  for  elliptic  interface  problems  with  piecewise  constant  coefficient  [6,  14],  and 
the  fast  algorithms  for  Poisson  and  biharmonic  equations  on  irregular  domains  [2,  7,  17,  16]. 
However,  the  augmented  approach  proposed  here  is  for  a  more  difficult  problem  and  it  may 
be  the  only  to  get  a  second  order  sharp  interface  method  to  solve  the  Stokes  flow  with 
discontinuous  viscosity.  The  key  to  the  success  of  an  augmented  approach  is  the  choice  of  the 
augmented  variable(s)  and  equation(s)  which  is  a  research  process  and  depends  on  problems. 

We  assume  that  the  domain  H  is  a  rectangle:  [a,  b]  x  [c,  d\ .  We  use  a  uniform  Cartesian 
grid 

Xi  =  a  +  ihx,  7  =  0,1,---  ,  m, 

Vj  =  c  +  jhy ,  J  =  0, 1,  •  •  •  ,n, 

so  that  a  fast  Poisson  solver  can  be  used  to  solve  the  three  equations  (3.17)-(3.19). 


b  —  a 
x  =  M  ’ 
d  —  c 

hy  =  ^T' 
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We  first  choose  a  set  of  control  points  {X/,.}  =  {(W,  I*,)},  k  =  1,  2,  •  •  •  ,  Nf,,  on  the  inter¬ 
face2.  The  auxiliary  variable  q(s)  =  (qi(s),  q2(s))  is  defined,  and  the  augmented  equations 
(3.20)  are  discretized,  at  {X*.}.  We  use  upper  case  letters  such  as  Pij,  Uij ,  Vtj,  Q^,,  for  the 
discrete  approximations  at  grid  points  and  at  those  control  points  on  the  interface,  respec¬ 
tively.  We  use  the  bold  face  upper  case  letters  without  subscripts  to  represent  the  vectors 
formed  by  those  discrete  components. 

Given  an  initial  guess  of  Q  at  the  control  points,  we  can  approximate  its  first  and  second 
order  tangential  derivatives  ^  and  ,  see  [11,  13].  Thus  all  the  jump  conditions  in  the 
system  (3.17)-(3.19)  are  known  and  we  can  solve  the  system  (3.17)-(3.19)  sequentially  using 
the  immersed  interface  method  [11,  12],  Note  that  the  jump  conditions  for  the  pressure  and 
the  velocity  now  are  decoupled  if  we  know  Q.  Since  the  solution  depends  on  Q,  the  solution 
can  be  written  as  P(Q),  U(Q)  and  V(Q). 

If  the  computed  U(Q)  and  V(Q)  satisfy  the  two  equations  in  (3.20),  then  P(Q), 
U(Q)/V  and  V(Q)//u  are  an  approximated  solution  to  the  original  system  (1.1)-(1.2). 
Otherwise,  we  use  a  third  order  accurate  linear  interpolation  scheme  to  evaluate  the  residual 
of  the  two  equations  in  (3.20)  which  will  be  explained  in  Section  3.2. 

The  core  of  our  algorithm  is  to  use  the  GMRES  iterative  method  [21]  to  solve  for  the 
augmented  variable  U  and  V  which  are  one  dimensional  vectors  defined  along  the  interface  T, 
or  more  precisely,  at  the  control  points  {X^,}.  However,  there  is  no  need  to  find  the  coefficient 
matrix  explicitly  as  we  will  explain  below.  Each  GMRES  iteration  requires  to  solve  three 
Poisson  equations  with  given  jumps  in  the  solution  and  its  normal  derivative,  and  a  least 
squares  interpolation  scheme  to  evaluate  the  residual.  It  is  worth  to  point  out  that  the  most 
expensive  part  of  the  algorithm,  which  is  to  solve  the  three  Poisson  equations,  can  be  done 
by  calling  a  fast  Poisson  solver  three  times. 

3.1  The  discrete  system  of  equations  in  the  matrix  vector  form 

Given  a  discrete  approximation  of  (<71,(72)  at  {X/J,  we  can  solve  the  first  three  equations 
(3.17)-(3.19)  using  the  immersed  interface  method  [12]  to  get  an  approximate  solution:  the 
pressure  P(Q),  the  scaled  velocity  U(Q)  and  V(Q).  Generally  the  computed  velocity 
(U(Q),  V(Q))  do  not  satisfy  the  two  augmented  equations  in  (3.20),  that  is,  (U,  V)  = 
(U//U,  V///)  may  not  be  continuous  across  the  interface. 

Let  us  put  the  discrete  solution  {Pij},  {Uij},  and  {14,-}  together  as  a  big  vector  VI  whose 
dimension  is  3 MN.  We  denote  also  the  vector  of  the  discrete  values  of  (q±,  (72)  at  the  control 
points  {Xfc}  by  Q  whose  dimension  is  2 N),.  Then  the  discrete  solution  of  (3.17)-(3.19)  given 

2In  a  front  tracking  method,  {Xr-}  is  a  set  of  points  that  represents  the  interface,  see  [13]  for  examples 
in  which  a  cubic  spline  is  used.  In  a  level  set  method,  the  control  points  can  be  taken  as  the  orthogonal 
projections  of  irregular  grid  points  on  the  interface,  see  [6,  14]  for  example. 
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Q  can  be  written  as 

AU  +BQ  =  F1  (3.21) 

for  some  vector  Fi  and  sparse  matrices  A  and  B.  It  requires  solving  three  Poisson  equations 
with  different  force  terms  and  jump  conditions  to  get  IA  . 

Once  we  know  the  solution  IA  given  Q,  we  can  use  (U,  V)  and  the  jump  conditions 
and  [^]  which  also  depend  on  Q,  to  get  [U(Q)]  =  [U(Q)//x]  and  [V(Q)]  =  [V(Q)//j] 
at  those  control  points  {Xfc},  1  <  k  <  Nf,.  If  both  ||  [U(Q)]  ||  and  |j  [V(Q)]  ||  are  smaller 
than  a  given  tolerance,  then  the  method  has  already  converged  and  Q,  U / //,  V//i  are  the 
approximate  solution.  The  interpolation  scheme  to  get  [U(Q)//u]  and  [V(Q)//r]  ,  which  will 
be  explained  in  detail  in  the  next  sub-section,  is  linearly  dependent  of  IA  ,  Q.  Therefore  we 
can  write 

[U(Q)]|  =  ([C//i]  ,  [u/m])T  =SU  +  £Q-F2,  (3.22) 

where  S  and  E  are  two  sparse  matrices,  and  F2  is  a  vector.  The  matrices  depend  on  the 
interpolation  scheme  and  are  only  used  for  theoretical  purpose  here  but  not  actually  con¬ 
structed  in  our  algorithm.  We  need  to  choose  such  a  vector  Q  that  the  continuity  condition 
for  the  velocity  is  satisfied  along  the  interface  T.  If  we  put  the  two  matrix- vector  equations 
(3.21)  and  (3.22)  together  we  get 


A  B  1  [  U 

S  E  Q 


(3.23) 


Note  that  Q  is  defined  only  on  a  set  of  points  {X^.}  on  the  interface  while  IA  is  defined  at 
all  grid  points.  The  Schur  complement  for  Q  is 

{E  -  SA~1B)Q  =  F2~  SA~1F1  =  F.  (3.24) 

If  we  can  solve  the  system  above  to  get  Q,  then  we  can  get  IA  easily.  Because  the  dimension 
of  Q  is  much  smaller  than  that  of  W  ,  we  expect  to  get  a  reasonably  fast  algorithm  for  the 
two-phase  Stokes  equation  if  we  can  solve  (3.24)  efficiently. 

In  implementation,  we  use  the  GMRES  [21]  to  solve  (3.24).  The  GMRES  method  only 
requires  the  matrix  vector  multiplication.  We  explain  below  how  to  evaluate  the  right  hand 
side  F  of  the  Schur  complement,  and  how  to  evaluate  the  matrix  vector  multiplication  needed 
by  the  GMRES  iteration.  We  can  see  why  we  do  not  need  to  form  the  coefficient  matrix 
E  —  SA~1B  explicitly. 


3.1.1  Evaluation  of  the  right  hand  side  of  the  Schur  complement 

First  we  set  Q  =  0  and  solve  the  de-coupled  system  (3.17)-(3.19),  or  (3.21)  in  the  discrete 
form,  to  get  IA  (0)  which  is  J4_1Fi  from  (3.21).  From  the  interpolation  scheme  (3.22),  we 
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also  have 

[U(0)]  =  SU{0)  +  E0-F2  =  SU  (0)-F2.  (3.25) 

Note  that  the  residual  of  the  Schur  complement  for  Q  =  0  is 

R{ 0)  =  (E  -  SA~1B)  0  —  F  =  — F 

=  -(F2-SA~1F1)  =  -F2  +  SU(0)  (3.26) 

=  [U(0)]|r 

which  gives  the  right  hand  side  of  the  Schur  complement  system  with  an  opposite  sign. 

3.1.2  Evaluation  of  the  matrix-vector  multiplication 

The  matrix-vector  multiplication  of  the  Schur  complement  system  given  Q  is  obtained  from 
the  following  two  steps: 

Step  1:  Solve  the  coupled  system  (3.17)-(3.19),  or  (3.21)  in  the  discrete  form,  to  get  tl  (Q). 

Step  2:  Interpolate  W  (Q)  using  (3.22)  to  get  [U(Q)]|r.  Then  the  matrix  vector  multipli¬ 
cation  is 

(E  -  SA~1B)Q=  [U(Q)]  —  [U(0)J  .  (3.27) 

This  is  because 

(E-SA~lB)  Q  =  EQ-SA^BQ 

=  EQ-S(^A~1F1  -W(Q))  (from  (3.21)  )  , 

=  EQ  +  S’W(Q)-F2  +  F2-^-1Fi 
=  [U(Q)]|  +F  ( from  (3.22)  )  , 

=  [U(Q)]|r  -  [U(0)]|r,  (  from  (3.26)  )  . 

Now  we  can  see  that  a  matrix  vector  multiplication  is  equivalent  to  solving  the  coupled 
system  (3.17)-(3.19),  or  (3.21)  in  the  discrete  form,  to  get  IA  ,  and  using  an  interpolation 
scheme  (3.22)  to  get  [U(Q)]|r  at  the  control  points. 

Since  we  know  the  right  hand  side  of  the  linear  system  of  equation  and  the  matrix- vector 
multiplication  of  the  coefficient  matrix,  it  is  straightforward  to  use  the  GMRES  or  other 
iterative  methods. 
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3.2  The  least  squares  interpolation  scheme  to  compute  the  residual. 

The  interpolation  scheme  (3.22)  to  evaluate  [U//i]  and  [V //j]  is  crucial  to  the  efficiency 
(accuracy  and  the  number  of  iterations  of  the  GMRES  iteration)  of  the  method.  To  reduce 
the  number  of  iterations,  it  is  important  to  couple  the  solutions  on  both  sides  of  the  interface 
using  the  jump  conditions.  Although  the  least  squares  interpolation  scheme  is  not  new  idea 
anymore,  the  details  vary  with  problems.  We  explain  the  least  squares  interpolation  scheme 
for  computing  (3.22)  to  see  why  we  have  the  second  matrix-vector  equation  in  expression 
(3.23). 

Given  an  approximation  to  the  augmented  variable  Q,  we  can  solve  for  the  pressure,  and 
then  the  velocity  (U,  V)  =  (U / g,  V//x)  from  (3.17)-(3.19).  Since 

ul  _  u+  u- 

M  J  h+ 

we  need  to  evaluate  {t/+}  and  {U~}  at  all  control  points  to  get  the  vector  [XJ/g].  To 
explain  the  idea,  however,  we  just  need  to  explain  the  interpolation  scheme  for  U~  (X)  at  a 
point  X  on  the  interface.  The  interpolation  scheme  can  be  written  as 

ks- 1 

U  (X)  =  ^2  Ik  Ui*+k,j*+k  +  C,  (3.28) 

k= o 

where  ks  is  the  number  of  grid  points  involved  in  the  interpolation  scheme,  ( Xi*,yj *)  is  the 
closest  grid  point  to  X,  and  C  is  a  correction  term.  We  should  point  it  out  that  a  one-sided 
interpolation  scheme  works  poorly  in  the  sense  that  the  convergence  speed  is  slow  for  the 
GMRES  iteration.  Below  we  discuss  how  to  determine  the  coefficients  and  the  correction 
term  C  using  the  information  from  both  sides  of  the  interface.  Note  that  jk  and  C  depend 
on  X.  But  for  simplicity  of  the  notation,  we  have  omitted  the  dependency. 

We  use  an  un-determined  coefficients  method  to  determine  the  coefficients  ^k  by  mini¬ 
mizing  the  truncation  error  of  (3.28)  when  J7j*+ k,j*+k  is  substituted  by  the  exact  solution 
u(xi*+k,  yj*+k).  Using  the  local  coordinates  system  centered  at  the  point  X,  see  (2.9),  and 
denoting  the  local  coordinates  of  (xi*+k,  Vj*+k)  as  (£&>%)>  we  have  the  following  from  the 
Taylor  expansion: 

u(xi*+k,  yj*+k)  =  u(£k,  Vk)  =  u±  +  Ckuf  +  r]ku± 

1  1  (3.29) 

+  2  +  ZkVkufv  +  ~  vltfv  + 

where  the  +  or  —  sign  is  chosen  depending  on  whether  (£&,%)  lies  on  the  +  or  —  side  of 
T,  u^,  u±,  •••,  u±v  are  evaluated  at  the  local  coordinates  (0,0),  or  X  =  (X,Y)  in  the 
original  coordinates  system.  Note  that  we  should  have  used  something  like  u(X,  Y )  =  u( 0, 0) 
to  distinguish  the  two  coordinate  systems.  However,  we  omit  the  bars  and  use  the  same 
notation  u(X,  Y)  =  u(0, 0)  for  simplicity. 
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We  carry  out  this  expansion  at  all  the  grid  points  used  in  the  interpolation  scheme  and 
plug  (3.29)  into  (3.28).  After  collecting  and  arranging  terms,  we  can  write 


u  “(X)  «  a\u  +  02  u+  +  d3  +  a4  u I"  +  a$  uv  +  a$ u+ 

+  O7  Ugg  +  08  +  Og  +  Oio  +  On  +  Oi2  U^C, 

where  the  Oj’s  are  given  by 


(3.30) 


Oi  = 

X 

7  k 

02  = 

2 

keK~ 

keK+ 

03  = 

X 

ikik 

04  = 

y,  ikik 

k&K~ 

keK+ 

0.5  = 

X 

Vklk 

06  = 

y  dkik 

keK~ 

keK+ 

«7  = 

IE 

ifok 

as  = 

2  &lk 

keK- 

- 

keK+ 

09  = 

^x 

vhk 

oio  = 

--  \  dhk 

k£K- 

- 

keK+ 

on  = 

=  X 

ikhklk 

012  = 

=  ^2  ^mik- 

k£l<- 

kei<+ 

+  _ 

U~  +  ( 

h  and  [f|] 

is  known  from 

(2.11).  From 

(3.31) 


following  interface  relations3: 


~+  —  .  dqi 

uv=uv+  — , 


dr\ 


^7777  -  ^7777  “I-  ^ 


du 

dn 

du 

dn 


d2(ll  r  , 

-w+lPx~9l] 


+ 


~+  —  ,  dqi  d 

ain  =  \n  +  ^  +  ^ 


d2qi 
dr 12  ’ 
du 
dn 


(3.32) 


where  k  is  the  curvature  of  the  interface  at  X.  Therefore  we  can  express  all  the  quantities 
from  +  side  in  (3.30)  in  terms  of  those  from  —  side  and  the  known  quantities.  Thus  (3.28), 

3Note  that  qi  now  is  a  quantity  defined  only  on  the  interface,  and  we  can  only  take  its  derivatives  along 
the  interface  which  is  called  surface  derivative  in  the  literature. 
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when  Ui*+k,j*+k  is  substituted  by  the  exact  solution  u(xi*+kiUj*+k )>  can  be  written  as 
i-(X)  «  £  Tfc  'U'{A'i*+ki  Vj*+k)  C 

k 

=  a\  u~  +  a2  u+  +  a3  +  a4  u ^  +  a5  u~  +  aeu+  +  07  u ^  +  a8 
+  ag  +  aio  Uqjj  +  an  +  ai2  u ^  +  C 

=  (ai  +  a2)tT  +  (a3  +  a4)u^  +  (a5  +  a6)u^  +  (07  +  a8)u^ 

+  (a9  +  aw)u~n  +  (an  +  ai2)u^  +  a2[u]  +  a4[%]  +  a6[u,?] 

+  a8[u^]  +  aiofu^]  +  ai2[u^,?]  +  C. 


To  minimize  the  local  truncation  error,  we  should  set  the  following  linear  system  of  equations 
for  the  coefficients  7*,  by  matching  the  terms  of  u~ ,  igf ,  •  •  • ,  u^,tl  : 


a\  +  a2  =  1, 

a3  +  a4  =  0, 

as  +  a6  =  0, 
07  +  a8  =  0, 
ag  +  aio  =  0, 

an  +  ai2  =  0. 


(3.33) 


The  system  of  equations  for  {7^.}  is  independent  of  jumps  which  means  we  can  calculate  {7*.} 
outside  of  the  GMRES  iteration.  Once  we  have  the  coefficients,  the  correction  term  is 


C  =  -  (  a2[u]  +  a4[%]  +  ae[uv]  +  a8[%]  +  ai0[a^]  +  ai2[u?,,]) 

d2[u\ 


du 

dqi  ( 

du 

a2  qi  —  a4 

dn 

dn 

k  — 


dr f 


+  [Px  -  gi  ] 


(3.34) 


aio 


(  d?q\ 

du 

\  ( <Au]  d 

du 

Cl. 
l3  , 
to 

1 

dn 

k  ai2  ,  K+  , 

J  \  dr]  dr] 

dn 

We  choose  a  neighborhood  of  {X}  that  contains  more  than  six  different  grid  points  so 
that  we  have  an  under-determined  system.  In  our  numerical  tests,  we  choose  ks  =  12,  that 
is,  we  selected  12  closest  grid  points  to  X  =  (X,  Y)  as  the  interpolation  stencil.  We  use 
the  singular  value  decomposition  (SVD)  to  find  the  least  squares  solution  which  also  has  the 
least  1-2  norm  among  all  the  solutions.  In  this  way,  the  magnitude  of  the  coefficients  7*, 
is  controlled  and  balanced.  The  least  squares  interpolation  plays  an  important  role  in  the 
stability  of  the  algorithm. 

The  only  trade-off  of  the  least  squares  interpolation  is  that  we  have  to  solve  an  under¬ 
determined  system  of  equations.  However,  the  size  of  the  linear  system  is  small  and  the 
coefficients  can  be  pre-determined  before  the  GMRES  iteration.  The  extra  time  needed  in 
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dealing  with  the  interface  is  usually  less  that  5%  of  the  total  CPU  time  and  the  percentage 
decreases  as  the  mesh  size  ( M  and  N)  increases. 

Remarks:  By  setting  a\  +  <12  =  0  and  0,3  +  <24  =  1  while  keeping  other  equations 
unchanged  in  (3.33),  we  can  easily  get  the  normal  derivative  of  the  solution  u.  This  is  the 
method  that  we  used  in  Section  4  for  the  grid  refinement  analysis  of  the  computed  normal 
derivative  of  the  velocity. 

There  are  also  a  few  issues  about  how  to  evaluate  px  and  py  at  grid  points  near  the 
interface.  We  also  need  to  evaluate  [px]  and  [py]  at  some  points  on  the  interface.  We  refer 
the  reader  to  [12]  for  the  detailed  information.  We  just  need  to  compute  those  quantities 
near  or  on  the  interface  to  first  order  accuracy  since  they  appear  at  the  right  hand  side  of 
the  Poisson  equations. 


4  Results  and  analysis  of  numerical  experiments 


We  first  present  several  grid  refinement  analysis  for  the  examples  in  which  we  know  the  exact 
solutions  to  check  the  order  of  accuracy.  Most  of  the  computations  are  done  on  workstations 
or  Laptop  PC’s  within  a  few  seconds  to  a  few  minutes  for  stationary  Stokes  problems.  It 
is  quite  challenging  to  construct  the  exact  solutions  for  incompressible  Stokes  flow  with  an 
interface.  Throughout  this  section,  the  computational  domain  is  P  =  [—2,  2]  x  [—2,  2]. 

Example  4.1 

We  start  with  a  simple  example  where  the  velocity  is  smooth  and  the  pressure  is  discontinuous 
across  the  interface.  The  exact  velocity  and  the  pressure  are  given  by 

u  =  y{x2  +  y2-  1)  ,  (x,y)eCl, 
v  =  -x  (x2  +  y2  -  l)  ,  (x,  y)  G  Cl, 

f  1,  if  x2  +  y2  <  1, 

p  =  < 

[0,  if  x2  +  y2  >  1. 

The  interface  is  the  unit  circle.  The  viscosity  is 

1,  if  x2  +  y2  <  1, 


p  =  <  1 

-,  if  x2  +  y2  >  1 . 

The  bounded  external  forcing  term  g  is  given  by 

f  -8 y,  if  x2  +  y2  <  1, 


9i  = 


92  = 


-4 y,  if  x2  T  y2  >  1, 

8x,  if  x2  +  y2  <  1, 
Ax,  if  x2  +  y2  >  1, 


(4.35) 

(4.36) 

(4.37) 


(4.38) 


(4.39) 


(4.40) 
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which  has  a  finite  jump  across  the  interface.  The  normal  and  tangential  force  density  are 


h 


h 


M-2 


<9u 


=  -l, 


<9u 

<9u 

»dn  '  T 

^  ' n 

=  -l, 


(4.41) 

(4.42) 


calculated  from  (2.4)  and  (2.6)  respectively. 

In  Table  1,  we  show  the  result  of  the  grid  refinement  analysis.  Since  the  exact  solution 
is  not  periodic,  we  use  the  Dirichlet  boundary  condition  when  we  solve  the  three  Poisson 
equations  (3.17)-(3.19)  so  that  we  can  check  the  accuracy  of  the  computed  solution.  This  did 
not  cause  any  inconsistence  since  we  know  that  the  solution  exists.  Similarly,  the  Dirichlet 
boundary  condition  is  used  for  other  examples  with  the  exact  solution  that  is  not  periodic  in 
this  section  based  on  the  same  arguments.  If  we  choose  0  as  the  initial  guess,  the  GMRES 
iteration  converges  quickly  since  0  is  the  exact  solution  of  the  augmented  variables.  Therefore 
we  choose  random  numbers  between  0  and  1  as  the  components  of  the  initial  guess  to  get  a 
realistic  evaluation  of  the  algorithm.  The  tolerance  for  the  GMRES  iteration  is  taken  as  10 ~6. 
The  errors  in  Table  1  are  measured  in  the  maximum  norm  at  all  grid  points,  for  example 


EU(N) 


1 

2 


max 

o  <i,j<N 


I  Uij 


—  u(xi,  Dj)\  +  max  \Vij  -  v(xi:  yj)\ 
0<i,j<N 


(4.43) 


where  u(xi,yj)  is  the  exact  solution  at  ( Xi,yj )  while  Uij  is  the  approximate  solution  and  so 
on.  In  all  the  tables  in  this  section,  N  is  the  number  of  grid  lines  in  both  x-  and  y-  directions. 
The  ratio 


p-order 


log  {Ep{N)/ EP(2N)) 
log  2 


(4.44) 


for  example,  is  an  indication  of  the  order  of  accuracy  for  the  pressure.  We  can  see  roughly 
second  order  convergence  for  all  the  quantities.  The  last  column  is  the  number  of  iterations 
(No.)  for  the  GMRES  method.  We  can  see  that  the  number  of  iterations  remains  roughly 
the  same  as  we  double  the  grid  lines  in  each  direction. 


Table  1:  Numerical  results  and  convergence  analysis  for  Example  4.1. 


N 

Ep 

p-order 

Eu 

u-order 

Edu/dn 

-order 

No. 

32 

8.2569  x  10"3 

6.5931  x  10"3 

2.4605  x  10“2 

10 

64 

3.0540  x  10"3 

1.4349 

1.7372  x  10"3 

1.9242 

6.0768  x  10"3 

2.0176 

11 

128 

9.4745  x  10-4 

1.6886 

3.9504  x  10-4 

2.1367 

1.6713  x  10-3 

1.8623 

12 

256 

2.6827  x  10-4 

1.8204 

8.2274  x  10-5 

2.2635 

4.1920  x  10-4 

1.9953 

13 

512 

7.4298  x  10"5 

1.8523 

2.5053  x  10"5 

1.7155 

1.0445  x  10"4 

2.0048 

14 
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Figure  2:  Linear  regression  analysis  for  the  pressure  and  the  velocity  for  Example  4.1.  The 
order  of  accuracy  for  the  pressure  and  the  velocity  are  1.9187  and  1.9416  respectively. 

Using  fixed  grids  to  solve  an  interface  problem,  the  errors  usually  do  not  decrease  rnono- 
tonically  as  we  refine  the  grid,  see  [14].  To  be  more  precise,  we  use  the  linear  regression 
analysis  to  find  approximate  order  of  accuracy.  In  Fig.  2,  we  show  the  error  plot  in  log-log 
scale  for  the  pressure  and  the  velocity  versus  the  grid  spacing  h  =  hx  =  hy.  which  shows 
that,  from  the  slopes,  the  average  convergence  order  of  the  pressure  and  the  velocity  are 
1.9187  and  1.9416,  respectively.  The  mesh  size  varies  from  N  =  100  to  N  =  320  according 
to  N  =  100  +  5k,  k  =  0,1,  •  •  •  ,44.  The  order  of  accuracy  for  the  normal  velocity  du/dn  is 
2.1928  from  the  linear  regression  analysis. 

Example  4.2 

In  the  first  example,  while  the  pressure  is  discontinuous,  the  velocity  is  smooth  and  vanishes 
at  the  interface.  Therefore  the  exact  solution  to  the  augmented  variable  q  =  [pu\  is  zero 
which  may  make  the  problem  easier  to  compute.  In  this  example,  we  keep  all  the  solution 
inside  the  unit  circle  unchanged,  but  set  both  pressure  and  the  velocity  outside  of  the  circle 
to  be  zero.  Therefore,  the  periodic  boundary  conditions  are  satisfied.  The  velocity  is  non¬ 
smooth  and  the  jump  in  the  normal  velocity  is  not  a  constant  along  the  interface.  Now  the 
normal  force  density  is  still  fi  =  —1,  but  the  tangential  force  density  is 

t  r  du  1  r  du  1 

h  =  -  ^  •  r  -  n—  •  n  =  -2.  (4.45) 

We  also  have  g  =  0  outside  of  the  circle  and  there  is  a  finite  jump  in  g  as  well.  In  Table  2,  we 
show  the  grid  refinement  analysis.  We  see  again  roughly  average  second  order  convergence 


.5  -3  -2.5  -2  -1.5  -1 

logloh 
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for  all  quantities.  The  number  of  iterations  is  small  and  does  not  change  much.  The  results 
of  the  linear  regression  analysis  for  the  pressure  and  the  velocity  are  given  in  Fig.  3  which 
confirms  average  second  order  accuracy  for  both  the  pressure  and  the  velocity. 


Table  2:  Numerical  results  and  convergence  analysis  for  Example  4.2. 


N 

Ep 

p-order 

Eu 

u-order 

Egu/dn 

Srorder 

No. 

32 

8.4430  x  10"3 

3.4549  x  10"3 

2.8308  x  10“2 

9 

64 

2.8405  x  10-3 

1.5716 

8.8800  x  10-4 

1.9600 

6.1798  x  10-3 

2.1956 

10 

128 

8.0952  x  10-4 

1.8110 

2.2666  x  10-4 

1.9700 

1.8260  x  10-3 

1.7589 

11 

256 

2.5417  x  10“4 

1.6713 

4.7693  x  10"5 

2.2487 

5.3612  x  10"4 

1.7681 

12 

512 

1.4086  x  10"5 

2.1296 

1.4086  x  10"5 

1.7595 

1.2538  x  10"4 

2.0962 

13 

Figure  3:  Linear  regression  analysis  for  Example  4.2.  The  average  convergence  order  for  the 
pressure  and  the  velocity  are  porder  =  1-9168,  uorder  =  2.1519. 


Example  4.3 

In  previous  examples,  the  force  density  are  constants,  and  [/i  ( u  •  n)]  =  0.  In  this  ex¬ 

ample,  we  construct  the  exact  solutions  in  such  a  way  that  all  the  jumps  and  their  derivatives 
along  the  interface  are  non-constant  functions.  The  exact  velocity  and  the  pressure  are  given 
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if  x2  +  y2  <  1, 
j  (®2  +  y2),  if  x2  +  y2  >  1, 


-^(l-x2),  if  x2  +  y2  <  1, 


if  x2  +  y2  >  1, 


3  o  3  \  ,  9  9 

--^cr  +  -x\y,  ifx2  +  y2<l, 


(4.46) 


(4.47) 


(4.48) 


if  x2  +  y2  >  1. 


The  bounded  external  forcing  term  g  is 


9i  = 


92  = 


(~Ix2  +  ^)y’  if  x'2  +  y2  < 

—2y+  y  if  x2  +  y2  >  1, 

3  o  3  3  LA  .  n  o  9 

_4X  +8X - 2~X’  lfx'  +y  <  ^ 

9  .  r  o  2  1 

—X,  if  x^  +  y  >  1, 


(4.49) 


(4.50) 


which  is  discontinuous  across  the  interface.  The  force  density  corresponding  to  the  singular 
Dirac  delta  function  in  the  normal  and  tangential  directions  are 


fi  =  ^  cos3  6  —  ^  cos  sin  6  —  ^  [y\  cos3  9  sin  8, 

h  =  \h+  +  |  [lA  cos2  0(1-2  cos2  9)  , 


(4.51) 

(4.52) 


respectively.  All  the  jump  conditions  (2.4)-(2.7)  are  satisfied  except  for  (2.5).  In  our  numerical 
test,  we  use  a  more  general  jump  condition 


—  =[g-n]  +  —  /2  +  2  7tTww(u  ■  n)  +  w 


(4.53) 


The  added  function  w  does  not  alter  the  difficulty  and  the  nature  of  the  problem  and  our 
algorithm  but  enables  us  to  check  the  order  of  accuracy  of  the  numerical  method.  In  this 
example,  w  is  given  by 


w  =  2  [//]  cos3  9  sin  9. 


(4.54) 


We  use  the  exact  Dirichlet  boundary  condition  for  the  pressure  and  the  velocity. 
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In  Table  3,  we  show  the  grid  refinement  analysis  for  different  jump  in  g.  We  scale  the 
problem  such  that  rna x{//“,/i+}  =  1.  We  test  our  results  for  / /i+  =  10,  10~3,  and  103. 
While  the  accuracy  does  depend  on  [/z],  the  average  convergence  rates  are  about  the  same 
(second  order  accurate).  Note  that,  there  are  two  very  different  scales  for  the  problems  in 
Table  3  (b)  and  Table  3  (c).  The  number  of  iterations  seems  to  be  dependent  of  the  jump  in 
l_i  but  not  the  mesh  size  N. 


Table  3:  Numerical  results  and  convergence  analysis  for  Example  4.3. 


(a)  n  =  1,  g+  =  0.1. 


N 

Ep 

J9-order 

Eu 

u-order 

Edu/dn 

gi-order 

No. 

32 

1.1252  x  10-2 

1.0331  x  10-2 

1.2397  x  10_1 

10 

64 

3.1248  x  10-3 

1.8484 

2.2832  x  10-3 

2.1779 

2.5154  x  10-2 

2.3011 

11 

128 

8.9338  x  10"4 

1.8064 

5.5784  x  10“4 

2.0331 

7.0234  x  10“3 

1.8405 

10 

256 

2.4296  x  10"4 

1.8786 

1.1291  x  10“4 

2.3047 

1.4007  x  10“3 

2.3260 

9 

512 

5.5515  x  10-5 

2.1297 

2.8135  x  10-5 

2.0047 

4.0508  x  10-4 

1.7899 

8 

(b)  fjL~  =  0.001,  //+  =  1. 


n 

Ep 

J9-order 

Eu 

u-order 

Edu/dn 

|))-order 

No. 

32 

1.4537  x  10-2 

9.9572  x  lO”1 

1.3266 

12 

64 

3.8694  x  10-3 

1.9095 

2.2471  x  10_1 

2.1477 

2.8381 

2.2247 

12 

128 

1.0974  x  10“3 

1.8180 

5.7087  x  10"2 

1.9768 

7.8753  x  lO”1 

1.8495 

11 

256 

3.4110  x  10~4 

1.6858 

1.1255  x  10"2 

2.3426 

1.8409  x  10“4 

2.0969 

11 

512 

6.6231  x  10-5 

2.3646 

2.8277  x  10-3 

1.9929 

4.7873  x  10-2 

1.9431 

9 

(c)  fi~  =  l,n+  =  0.001. 


n 

Ep 

J9-order 

Eu 

u-order 

Edu/dn 

an'Order 

No. 

32 

2.8939  x  10-2 

2.2153 

1.7133 

27 

64 

5.4693  x  10-3 

2.4036 

3.6865  x  lO”1 

2.5872 

2.5522  x  10_1 

2.7470 

25 

128 

1.6101  x  10“3 

1.7642 

9.1288  x  10"2 

2.0138 

5.7709  x  10“2 

2.1449 

25 

256 

4.1020  x  10~4 

1.9728 

2.4121  x  10“2 

1.9201 

1.4379  x  10“2 

2.0048 

31 

512 

1.3470  x  10~4 

1.6066 

6.7498  x  10-3 

1.8374 

3.3450  x  10-3 

2.1039 

30 

Example  4.4 


As  a  final  test,  we  present  an  example  of  moving  interface  problem.  The  driven  force  is  the 
surface  tension,  that  is,  the  normal  force  is  given  by  fi  =  k,  where  n  is  the  curvature  of  the 
interface  T.  The  tangential  force  is  zero.  The  interface  T  moves  with  the  same  velocity  as  the 
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fluid  surrounding  it,  i.e.,  -jj-  =  u.  The  set-up  is  almost  the  same  as  the  problem  described  in 
[12]  except  that  the  viscosity  is  now  discontinuous.  Since  the  emphasis  of  this  paper  is  about 
the  new  method  for  stationary  Stokes  equations,  we  omit  the  details  about  the  algorithm  for 
the  moving  interface  problems,  but  just  show  one  simulation  result.  The  initial  interface  in 
polar  coordinates  is  set  to  be 

p  =  1  +  0.2 sin(5a),  0  <  a  <  2ir.  (4.55) 

The  interface  will  move  to  its  equilibrium,  a  circle.  In  Fig.  4,  we  show  the  initial  interface 
and  a  computed  approximation  to  the  equilibrium.  We  tested  three  different  cases  with 
(p~ ,  p+)  =  (1,0.01),  (1,1),  and  (0.01,1).  Fig.  5  shows  the  close-up  plot  of  the  computed 
interface  at  t  =  0.5  for  the  three  different  cases.  As  it  is  expected  that  the  interface  moves 
faster  if  the  viscosity  is  small.  How  the  jump  in  the  viscosity  affects  the  motion  of  the  interface 
is  still  under  investigation. 


2 
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1 
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Figure  4:  Plots  of  the  initial  interface  (the  dotted  line)  and  the  final  equilibrium  (the  solid 
line)  of  Example  4.4. 
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Conclusion  and  acknowledgments 

In  this  paper,  a  new  second  order  accurate  numerical  method  has  been  developed  for  in¬ 
compressible  stationary  Stokes  equations  with  discontinuous  viscosity  in  which  the  jump 
conditions  for  the  pressure  and  the  velocity  are  coupled  together.  The  idea  is  to  introduce 
two  augmented  variables  that  are  only  defined  along  the  interface  so  that  the  jump  conditions 
can  be  decoupled.  The  GMRES  iterative  method  then  is  used  to  solve  the  Schur  complement 


Augmented  techniques  for  Stokes  flows 


21 


(a). 


(b). 


Figure  5:  Plots  of  the  computed  interface  near  the  two  points  A  and  B  with  different  viscosity 
ratio  at  t  =  0.2.  The  dotted  line  is  the  initial  interface.  The  dash-dotted  line  is  the  result 
with  (g- ,  fi+)  =  (1,  1).  The  solid  line  is  the  result  with  (/x“,  /i+)  =  (0.01,  1).  The  dashed 
line  is  the  result  with  ,  /j,+)  =  (1,  0.01). 

system  for  the  augmented  variables.  The  main  cost  in  one  step  GMRES  iteration  is  about 
three  calls  to  a  fast  Poisson  solver.  Numerical  examples  shown  the  efficiency  of  the  method  in 
accuracy  and  the  speed.  We  believe  that  the  new  method  may  be  the  first  second  order  sharp 
interface  method  for  Stokes  flow  with  discontinuous  viscosity.  The  idea  should  be  applicable 
to  other  interface  problems  with  coupled  jump  conditions. 
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